if (!require("pacman"))
install.packages("pacman")
# Use pacman::p_load to install and load CRAN packages
pacman::p_load(
countdown,
# datasauRus,
dplyr,
elevatr,
gganimate,
ggplot2,
ggthemes,
ggspatial,
gifski,
glue,
grid,
gridExtra,
here,
hms,
knitr,
lubridate,
metR,
openintro,
readr,
rnaturalearth,
rnaturalearthdata,
rnaturalearthhires,
terra,
scales,
sf,
tidyverse,
viridis
)
# Handle GitHub package separately
if (!require("dsbox")) {
# Install devtools if not present
if (!require("devtools"))
install.packages("devtools")
devtools::install_github("tidyverse/dsbox")
library(dsbox)
}
ggplot2::theme_set(ggplot2::theme_light(base_size = 12))
options(width = 65)
knitr::opts_chunk$set(
fig.width = 7, # 7" width
fig.asp = 0.618, # the golden ratio
fig.retina = 3, # dpi multiplier for displaying HTML output on retina
fig.align = "center", # center align figures
dpi = 300 # higher dpi, sharper image
)
save_figs <- TRUEProject Question 2
0 - Setup
vesuvius_data <- read.csv(here("data","vesuvius.csv"))
vesuvius_data_filtered <- vesuvius_data |>
drop_na(
latitude,
longitude,
depth_km,
duration_magnitude_md
) |>
mutate(
depth_boundary_1km = case_when(
depth_km < 1 ~ "< 1 km",
depth_km < 2 ~ "< 2 km",
depth_km < 3 ~ "< 3 km",
depth_km < 4 ~ "< 4 km",
depth_km < 5 ~ "< 5 km",
TRUE ~ "> 5 km"),
depth_boundary_1km = factor(
depth_boundary_1km,
levels = c("< 1 km",
"< 2 km",
"< 3 km",
"< 4 km",
"< 5 km",
"> 5 km")),
depth_boundary_0.5km = case_when(
depth_km < 0.5 ~ "< 0.5 km",
depth_km < 1 ~ "0.5 ~ 1 km",
depth_km < 1.5 ~ "1 ~ 1.5 km",
depth_km < 2 ~ "1.5 ~ 2 km",
depth_km < 2.5 ~ "2 ~ 2.5 km",
depth_km < 3 ~ "2.5 ~ 3 km",
depth_km < 3.5 ~ "3 ~ 3.5 km",
depth_km < 4 ~ "3.5 ~ 4 km",
depth_km < 4.5 ~ "4 ~ 4.5 km",
depth_km < 5 ~ "4.5 ~ 5 km",
TRUE ~ "> 5 km"),
depth_boundary_0.5km = factor(
depth_boundary_0.5km,
levels = rev(c("< 0.5 km",
"0.5 ~ 1 km",
"1 ~ 1.5 km",
"1.5 ~ 2 km",
"2 ~ 2.5 km",
"2.5 ~ 3 km",
"3 ~ 3.5 km",
"3.5 ~ 4 km",
"4 ~ 4.5 km",
"4.5 ~ 5 km",
"> 5 km")))
) |>
filter(
!is.na(depth_boundary_0.5km)
) |>
mutate(
depth_rescaled = rescale(depth_km, to = c(1, 11)),
mag_rescaled = rescale(duration_magnitude_md, to = c(0.4, 0.9))
) |>
arrange(
depth_boundary_0.5km,
desc(depth_rescaled)
) |> glimpse()Rows: 8,475
Columns: 15
$ event_id <int> 34136, 5595, 13138, 1499, 21103, …
$ time <chr> "2024-04-14T16:05:54Z", "2021-08-…
$ latitude <dbl> 40.86050, 40.82580, 40.82570, 40.…
$ longitude <dbl> 14.34917, 14.43120, 14.43000, 14.…
$ depth_km <dbl> 9.35, 5.11, 4.78, 4.42, 4.00, 4.0…
$ duration_magnitude_md <dbl> 2.30, 1.60, 2.18, 2.50, 1.70, 0.0…
$ md_error <dbl> 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3…
$ area <chr> "Mount Vesuvius", "Mount Vesuvius…
$ type <chr> "earthquake", "earthquake", "eart…
$ review_level <chr> "preliminary", "preliminary", "pr…
$ year <int> 2024, 2021, 2021, 2019, 2014, 201…
$ depth_boundary_1km <fct> > 5 km, > 5 km, < 5 km, < 5 km, <…
$ depth_boundary_0.5km <fct> > 5 km, > 5 km, 4.5 ~ 5 km, 4 ~ 4…
$ depth_rescaled <dbl> 11.000000, 6.460385, 6.107066, 5.…
$ mag_rescaled <dbl> 0.8069767, 0.7255814, 0.7930233, …max_depth <- max(vesuvius_data$depth_km, na.rm = TRUE)
min_depth <- min(vesuvius_data$depth_km, na.rm = TRUE)
print(paste("Maximum depth (km):", max_depth))[1] "Maximum depth (km): 9.35"print(paste("Minimum depth (km):", min_depth))[1] "Minimum depth (km): 0.01"max_mag <- max(vesuvius_data$duration_magnitude_md, na.rm = TRUE)
min_mag <- min(vesuvius_data$duration_magnitude_md, na.rm = TRUE)
print(paste("Maximum Magnitude:", max_mag))[1] "Maximum Magnitude: 3.1"print(paste("Minimum Magnitude:", min_mag))[1] "Minimum Magnitude: -2"original <-
ggplot(
vesuvius_data,
aes(x = longitude,
y = latitude)
) +
geom_point(
aes(color = duration_magnitude_md,
alpha = abs(depth_km)),
size = 3
) +
scale_color_viridis_c(
option = "inferno",
name = "Magnitude"
) +
scale_alpha_continuous(
name = "Depth (km)",
range = c(0.1, 0.8)
) +
labs(
title = "Seismic Events at Mount Vesuvius",
subtitle = "Point Color by Magnitude\nPoint Transparency by Depth",
x = "Longitude",
y = "Latitude",
caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
)
original
if (save_figs) {
ggsave("images/orignal_vesuvius_plot.png", original)
}histogram <-
ggplot(
vesuvius_data_filtered,
aes(x = duration_magnitude_md)
) +
geom_histogram(
binwidth = 0.2,
fill = 'red',
color = "white"
) +
facet_wrap(
~ depth_boundary_0.5km,
scales = "free_y"
) +
labs(
title = "Histogram of Duration Magnitude",
subtitle = "Distribution Facets by Depth",
x = "Magnitude (md)",
y = "Count",
caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
) +
theme(
plot.title.position = "plot"
)
histogram
if (save_figs) {
ggsave("images/histogram_vesuvius_plot.png", histogram)
}# Define Vesuvius bounding box (in WGS84 lon/lat)
vesuvius_bbox <- matrix(c(14.33, 40.78, 14.48, 40.89), ncol = 2, byrow = TRUE)
# Convert bbox to sf polygon
vesuvius_sf <- st_as_sf(
st_sfc(st_polygon(list(rbind(
c(vesuvius_bbox[1,1], vesuvius_bbox[1,2]),
c(vesuvius_bbox[2,1], vesuvius_bbox[1,2]),
c(vesuvius_bbox[2,1], vesuvius_bbox[2,2]),
c(vesuvius_bbox[1,1], vesuvius_bbox[2,2]),
c(vesuvius_bbox[1,1], vesuvius_bbox[1,2])
))), crs = 4326))
# Get elevation data from AWS Terrain Tiles (zoom 10)
elev_raster <- get_elev_raster(
locations = vesuvius_sf,
z = 10,
clip = "locations")
# Convert to dataframe for ggplot
elev_df <- as.data.frame(
elev_raster,
xy = TRUE)
colnames(elev_df) <- c("x", "y", "elevation")static_plot <- ggplot() +
metR::geom_contour_fill(
data = elev_df,
aes(x = x,
y = y,
z = elevation/1000,
group = NULL),
bins = 10
) +
scale_fill_gradientn(
colors = gray.colors(10, start = 1.0, end = 0.3),
name = "Elevation (km)",
limits = c(-0.111, 1.246), # Match elevation range from elev_df (km)
breaks = seq(0, 1.2, by = 0.3)
) +
geom_point(
data = vesuvius_data_filtered,
aes(x = longitude,
y = latitude,
color = duration_magnitude_md,
size = depth_rescaled),
alpha = 0.75
) +
scale_size_continuous(
name = "Depth (km)",
breaks = c(2.5, 5, 7.5),
range = c(1, 11) # adjust to control point size
) +
scale_color_viridis_c(
option = "magma",
name = "Magnitude (md)"
) +
scale_x_continuous(
limits = c(14.34, 14.48),
expand = c(0, 0)
) +
scale_y_continuous(
limits = c(40.78, 40.88),
expand = c(0, 0)
) +
labs(
title = "Seismic Events at Mount Vesuvius",
subtitle = "Point Color by Magnitude\nPoint Size by Depth",
x = "Longitude",
y = "Latitude",
caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
) +
theme(
plot.title.position = "plot",
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
legend.direction = "horizontal",
legend.box = "vertical", # stack legends on top of each other
legend.box.just = "left",
legend.spacing.x = unit(0.3, "cm"), # tighter spacing
legend.spacing.y = unit(0.3, "cm"), # tighter spacing
legend.title = element_text(size = 12),
legend.text = element_text(size = 10)
) +
guides(
fill = guide_colorbar( # Elevation
direction = "horizontal",
title.position = "top",
barwidth = unit(2.8, "cm"),
barheight = unit(0.4, "cm"),
order = 1),
color = guide_colorbar( # Magnitude
direction = "horizontal",
title.position = "top",
barwidth = unit(2.8, "cm"),
barheight = unit(0.4, "cm"),
order = 2),
size = guide_legend( # Depth
direction = "horizontal",
title.position = "top",
label.position = "bottom",
override.aes = list(alpha = 0.7),
keywidth = unit(0.15, "cm"),
keyheight = unit(0.4, "cm"),
order = 3),
alpha = "none"
)
static_plot
if (save_figs) {
ggsave("images/static_vesuvius_plot.png", static_plot)
}ggplot2::theme_set(ggplot2::theme_light(base_size = 14))
options(width = 65)
knitr::opts_chunk$set(
fig.width = 7, # 7" width
fig.asp = 0.618, # the golden ratio
fig.retina = 2, # dpi multiplier for displaying HTML output on retina
fig.align = "center", # center align figures
dpi = 150 # higher dpi, sharper image
)vesuvius_anim_plot <- ggplot() +
geom_point(
data = vesuvius_data_filtered,
aes(x = longitude,
y = latitude,
color = duration_magnitude_md,
size = depth_rescaled),
alpha = 0.75
) +
scale_x_continuous(
limits = c(14.34, 14.48),
expand = c(0, 0)
) +
scale_y_continuous(
limits = c(40.78, 40.88),
expand = c(0, 0)
) +
scale_size_continuous(
name = "Depth (km)",
breaks = c(2.5, 5, 7.5),
range = c(1, 11)
) +
scale_color_viridis_c(
option = "magma",
name = "Magnitude (md)"
) +
labs(
title = "Seismic Events at Mount Vesuvius",
subtitle = "Point Color by Magnitude\nPoint Size by Depth: {closest_state}",
x = "Longitude",
y = "Latitude",
caption = paste(
"Source: Seismic data for Mount Vesuvius",
"TidyTuesday 2025-05-13",
"Graphic: Weston Scott",
sep = "\n"
)
# caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
) +
theme(
plot.title.position = "plot",
plot.title = element_text(margin = margin(b = 23),
size = 18),
axis.text = element_text(size = 12),
legend.direction = "horizontal",
legend.box = "vertical", # stack legends on top of each other
legend.box.just = "left",
legend.spacing.x = unit(0.3, "cm"), # tighter spacing
legend.spacing.y = unit(0.3, "cm"), # tighter spacing
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
plot.caption = element_text(margin = margin(t = 30), size = 10, hjust = 1)
) +
guides(
color = guide_colorbar(
direction = "horizontal",
title.position = "top",
barwidth = unit(2.8, "cm"),
barheight = unit(0.4, "cm"),
order = 1),
size = guide_legend(
direction = "horizontal",
title.position = "top",
label.position = "bottom",
override.aes = list(alpha = 0.7),
keywidth = unit(0.15, "cm"),
keyheight = unit(0.4, "cm"),
order = 2),
alpha = "none"
) +
transition_states(
depth_boundary_0.5km,
transition_length = 3,
state_length = 2
) +
shadow_mark(
past = TRUE,
future = FALSE,
alpha = 0.2
) +
enter_fade() +
exit_fade()
gganimate::animate(vesuvius_anim_plot)
if (save_figs) {
gganimate::anim_save("images/vesuvius_animation.gif", vesuvius_anim_plot)
}